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We study cold atomic gases in a disorder potential and analyze the correlations between different 
systems subjected to the same disorder landscape. Such independent copies with the same disorder 
landscape are known as replicas. While in general these are not accessible experimentally in con- 
densed matter systems, they can be realized using standard tools for controlling cold atomic gases 
in an optical lattice. Of special interest is the overlap function which represents a natural order pa- 
rameter for disordered systems and is a correlation function between the atoms of two independent 
replicas with the same disorder. We demonstrate an efficient measurement scheme for the determi- 
nation of this disorder-induced correlation function. As an application, we focus on the disordered 
Bose-Hubbard model and determine the overlap function within perturbation theory and a numer- 
ical analysis. We find that the measurement of the overlap function allows for the identification of 
the Bose-glass phase in certain parameter regimes. 



PACS numbers: 03.75 Lm, 64.60 Cn, 42.50 -p 
I. INTRODUCTION 

The interplay between disorder and interaction gives 
rise to a plethora of fundamental phenomena in con- 
densed matter systems. The most predominant exam- 
ples include spin glasses [U El [3] , the superconducting-to- 
insulator transition in thin superconducting films [H [S] , 
and localization phenomena in fermionic systems such as 
weak localization and the metal-insulator transition [S]. 
While the nature of order in spin glasses and its theoreti- 
cal description is still highly debated [3 El [i H [101 [H [H 
[T51 [H] . a substantial contribution to the understanding 
of bosons in a disordered medium has been provided by 
work on the disordered Bose-Hubbard model introduced 
by Fisher et al. [15 . The disordered Bose-Hubbard model 
has recently attracted considerable attention due to its 
potential realization with cold atomic gases in an optical 

lattice [n [T71 nHiiiniiiniini mi- 

A general challenge in the description of glass phases 
in disordered media is the absence of a simple order 
parameter distinguishing the different ground states. 
This problem becomes evident in the disordered Bose- 
Hubbard model where the phase diagram is determined 
by the competition between a superfluid phase, the Mott- 
insulator, and a Bose-glass phase. While the super- 
fluid phase exhibits a finite condensate fraction and is 
characterized by off-diagonal long-range order, the Bose 
glass is only distinguished from the incompressible Mott- 
insulator by a vanishing excitation gap and a finite com- 
pressibility. However, in experiments on cold atomic 
gases the excitation gap and the compressibility are dif- 
ficult to determine accurately and are also obscured 
by the finite harmonic trapping potential. The chal- 
lenge is therefore to develop measurement schemes al- 




FIG. 1: (Colour online) Schematic representation of the im- 
plementation of disorder replicas with the cold atomic gases 
toolbox: Two planes a and /3 with equal disorder realisations, 
illustrated here for the case of a disorder landscape intro- 
duced using a second particle species (black dots), in which 
additional probe atoms (blue/lighter dots) evolve. The planes 
can be combined to measure correlation functions, such as the 
Edwards- Anderson overlap function. 



lowing for the experimental determination of these ob- 
servables or to develop new observables to characterize 
the glass phase. While inaccessible in "real materials," 
the Edwards- Anderson order parameter |T1 123] is com- 
monly studied analytically and measured numerically to 
quantify the "order" in a disordered system. It can be ex- 
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pressed as the correlation between independently evolv- 
ing systems with the same disorder landscape (so-called 
replicas) , or as a correlation of the same system at differ- 
ent temporal measurements. Because the latter depends 
on an extra variable (temporal measurement window) it 
is advantageous to study two replicas of the system with 
the same disorder. Thus, the measurement of this or- 
der parameter requires the preparation of several sam- 
ples with exactly the same disorder landscape, and the 
subsequent measurement of correlations between these. 

We demonstrate that such a procedure is feasible in 
cold atomic gases allowing one to gain access to char- 
acteristic properties of disordered systems naturally hid- 
den in condensed matter systems. The basic idea is to 
focus on cold atomic gases in an optical lattice: along 
one direction, the optical lattice is very strong and di- 
vides the system into independent two-dimensional (2D) 
planes [Ml EIJ [26l [27j. In addition, the system is sub- 
jected to a disorder potential and we are interested in 
the situation where in each plane the same disorder land- 
scape is realized (see figure [ij . Although the atoms in 
different planes are decoupled, the presence of the same 
disorder landscape within each plane induces a correla- 
tion between the different realizations: this correlation is 
of special interest in the replica theory of spin glasses and 
allows one to measure the overlap function, a character- 
istic property of spin glasses. Furthermore, we show that 
this correlation function is accessible in experiments on 
cold gases: The main idea is to quench the motion of the 
atoms by a strong optical lattice, and combine the differ- 
ent planes into a single one. Subsequently, the particle 
number occupation within each well carries the informa- 
tion about the correlation function. The possibility for 
the accurate determination of the particle number within 
each well of an optical lattice has recently been demon- 
strated experimentally for the superfluid-Mott-insulator 
quantum phase transition |28j . 

Note that this measurement scheme for the overlap 
between system with the same disorder can be applied 
to any disordered one- or two-dimensional system real- 
ized with cold atomic gases in an optical lattice. Of spe- 
cial interest is the realization of spin glasses and their 
study in the quantum regime. As an application, we 
focus here on the disordered Bose-Hubbard model and 
calculate the overlap function analytically in different 
regimes and compare it with one-dimensional numerical 
simulations. We find that the overlap function q exhibits 
a sharp crossover from the Mott-insulating phase with 
g « [25] to the Bose-glass phase with q « p{l — p), 
with p G (0,1), and thus makes it possible to distin- 
guish the two phases. Because the superfluid phase can 
be detected via the interference peaks in a time-of-flight 
measurement, we propose that this novel measurement 
scheme for the overlap function can be used for a quali- 
tative experimental verification of the phase diagram for 
the disordered Bose-Hubbard model. 

Different implementations of disorder or quasi-disorder 
in cold atomic gases have been discussed and experi- 



mentally realized. Several groups attempted to search 
for traces of Anderson localization and the interplay of 
disorder-nonlinear interactions in Bose-Einstein conden- 
sates (BECs). The first experiments [501 El ISH [Ml 131] 
were performed with laser speckles that had a disor- 
der correlation length larger than the condensate heal- 
ing length. Localisation effects were thus washed out by 
interactions. As an alternative, superlattice techniques — 
which are combinations of several optical potentials with 
incommensurable spatial periods — were used to produce 
quasi-disordered potentials with short correlations. This 
approach allowed the observation of some signatures of 
the Bose glass by the Florence group [3S]. Only very re- 
cently the Palaiseau group has managed to create speckle 
potentials on the sub-micron scale and directly observe 
Anderson localisations effects in a BEG released in a one- 
dimensional waveguide |36j based on the theoretical pre- 
dictions of [37l l38] . The Florence group reported ob- 
servations of localisations phenomena in quasi-disordered 
potentials in BEG of K'^^ , which allows for complete con- 
trol of the strength atomic interactions using Feshbach 
resonances |39j . It is worth noting that experimental at- 
tempts for local addressability in an optical lattice also 
offer the possibility to create disorder with correlation 
lengths comparable to the lattice spacing [40l |4T| . Mix- 
tures of cold gases, on the other hand, provide an alterna- 
tive approach for the realization of disorder by quench- 
ing the motion of one species by a strong optical lat- 
tice, which then provides local impurities for the second 
species [42, 43 . While production of 2D planes with 
equal disorder realisations follows naturally if the disor- 
der is implemented with laser speckles, we also show how 
such a system can be realized in the case of disorder in- 
duced by a mixture of atomic gases. 

In section [nj we start with a detailed description of the 
disordered Bose-Hubbard model and the different disor- 
der realizations. After a short review of the phase dia- 
gram of the disordered Bose-Hubbard model, we provide 
the definition of the overlap function q characterizing the 
disorder-induced correlations between independent real- 
izations of the system. In section jlllj we describe in de- 
tail the preparation of the system and the subsequent 
measurement scheme for the overlap function. In sec- 
tion IV we determine the overlap for the disordered Bose- 



Hubbard model via perturbation theory for physically- 
relevant limits and compare the result with numerical 
simulations of the one-dimensional system. Details of 
the calculations are presented in the appendices. 



II. BOSE-HUBBARD MODEL WITH 
DISORDER 



Gold atomic gases subjected to an optical lattice are 
well described by the Bose-Hubbard model [34^ ^ 

(ij) i i 
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with b] (bi) the bosonic creation (annihilation) operator 
and Tii = b\bi the particle number operator at the lattice 
site i. The first term describes the kinetic energy of the 
atoms with hopping energy J between nearest-neighbour 
sites, the second term accounts for the onsite interac- 
tion between the atoms, and the third term describes 
the disorder potential with random-site off-sets of the 
chemical potential ^. 

The disorder potential {Aj} can be generated via dif- 
ferent methods [Ml ISU [33 |42l |43] and is determined 
by a probability distribution P{5) describing the proba- 
bility of having an onsite shift with strength 5. The mean 
square of the disorder distribution 



d5 S^P{S) 



(2) 



gives rise to a characteristic energy scale A of the dis- 
order potential (note that the mean energy shift of the 
disorder potential is absorbed in the definition of the 
chemical potential). In addition, the disorder potential 
is characterized by a correlation length. Here, we are in- 
terested in short-range disorder with the disorder in dif- 
ferent wells independent of each other. The probability 
distribution P{S) and its correlation length depend on the 
source of the disorder potential, which varies depending 
on the microscopic implementation in cold gases. The 
most promising possibilities in order to produce disor- 
der with short correlation length involve the use of laser 
speckle patterns and mixtures of different cold atomic 
gases. While first experiments with laser speckles had a 
disorder correlation length larger than the lattice spacing 
[5ni 1^ . experimental efforts towards local addressabil- 
ity in an optical lattice also offer the possibility for the 
creation of disorder with a correlation length compara- 
ble to the lattice spacing, as well as Gaussian probability 
distributions [40l[41]. Alternatively, a disorder potential 
can also be created in mixtures of cold atomic gases in 
optical lattices [55', 15^ by quenching the motion of one 
species. Then, the disorder correlation length is of the 
order of the Wannier function, which can be smaller than 
the lattice spacing due to the strong lattice that quenches 
the motion of the disorder species, whilst the probability 
distribution becomes bimodal with A^ = ±A; here we 
are primarily interested in such a short ranged disorder 
with Ai being independent in the different wells. Note 
that both Gaussian disorder, as well as bimodal disor- 
der generally give rise to different physical phenomena in 
glass physics and thus both are interesting in their own 
right. 

The zero-temperature phase diagram of the Hamilto- 
nian (jlj was first studied by Fisher et al. where 
three different phases were discussed: the superfluid, the 
Mott-insulator, and the Bose-glass phases. In the two- 
dimensional regime of interest here, the superfluid phase 
appears for large hopping energies J > J7, A and is char- 
acterized by off-diagonal (quasi) long-range order (finite 
superfiuid condensate) and a linear excitation spectrum 
giving rise to a finite compressibility. On the other hand. 



for dominating interaction energies U^A,J the ground 
state corresponds to an incompressible Mott-insulator 
phase with an excitation gap and a commensurate fill- 
ing factor, i.e., the averaged particle number in a single 
well {rii) G N is an integer. While the quantum phase 
transition from the superfluid to the Mott-insulator has 
been experimentally identified f47], the disorder potential 
gives rise to an additional Bose-glass phase character- 
ized by a vanishing excitation gap and finite compress- 
ibility; see figure ^ for a sketch of the phase diagram. 
The details of the phase diagram have been studied via 
analytical methods in the regime of Anderson localiza- 
tion [331 HHl (SOI and in other regimes via numeri- 
cal methods such as quantum Monte Carlo [51], density 
matrix renormalization group (DMRG) [55], dynamical 
mean-field theory [53^ , and analytical methods [331 ■ 
Furthermore, recent interests focused also on the appear- 
ance of alternative phases for different disorder types in 
the Bose- Hubbard model, e.g., off-diagonal disorder giv- 
ing rise to a Mott-glass phase [55] while random onsite 
interactions can give rise to a Lifshits glass [57] . 
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FIG. 2: Sketch of the phase diagram for the disordered 
Bose-Hubbard model: the Mott-insulator (MI) appears for 
integer filling and the overlap function q vanishes for small 
hopping J/U ^ 1, while in the Bose-glass phase (BG) the 
overlap function approaches a finite value q = p{l — p) with 
p G (0, 1) characterizing the bimodal disorder distribution 
(see section IV I. Consequently a measurement of the overlap 



allows for a clear distinction between the Bose-glass and the 
Mott-insulator phase in this regime. In turn, the superfiuid 
phase (SF) is characterized by off-diagonal long-range order 
resulting in coherence peaks in a time of flight experiment. 

The superfluid phase in the absence of disorder is ex- 
perimentally identified via a measurement of the coher- 
ence peaks characterizing the condensate fraction, while 
the transition to a Mott-insulating phase is character- 
ized by the disappearance of these interference patterns 
and as well as a change in the behaviour of excitations 
[Ml H7] . However, the identification of the Bose-glass 
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phase in the presence of disorder requires an additional 
observable allowing for the distinction between the Mott- 
insulator and the Bose-glass, where the condensate frac- 
tion vanishes. Such an additional property is known from 
spin-glass theory as the Edwards- Anderson order param- 
eter jl] |23] which appears as an order parameter in the 
mean-field theory on spin glasses [58' . In the present sit- 
uation with a disordered Bose-Hubbard model, its gen- 
eralization leads to 

qEA = - nf]a,, = [{ni)\y - {[{ni)]avf , (3) 

where n = [(ri)^]^^ denotes the mean particle density in 
the sample. It is important to note the different aver- 
ages involved in the definition of the order parameter 
qsA'- for a fixed disorder realization, the average (...) de- 
notes the thermodynamic average over the ground state 
of the system, while [..]av describes the disorder average 
over different disorder realizations. While the experimen- 
tal determination of this quantity in the Bose-Hubbard 
model requires an exact state tomography of the system 
for each lattice site, a simpler and alternative route is ob- 
tained by preparing two systems with the same disorder 
landscape and measuring the correlations between these 
two decoupled systems: the remaining correlation is dis- 
order induced and disappears for weak disorder. This 
so-called overlap function q between the two systems re- 
alized with the same disorder landscape takes the form 

'Z-[«)("f)]a.-[«)]a.[(nf)]a. (4) 

and where a and [3 describe the two different systems 
with the same disorder. This definition follows in close 
analogy of the mean-field order parameter in the replica 
theory of spin glasses [58] . 

III. MEASUREMENT OF THE OVERLAP 
FUNCTION 

In this section we outline the key process of this work, 
i.e., the basic experimental procedure by which physical 
replicas can be prepared and used to measure the overlap 
function q defined in ^ . This consists of three essential 
steps: (1) initial preparation of disorder replicas, (2) in- 
troduction of probe atoms and (3) the measurement pro- 
cess, involving recombination of the replicas and spec- 
troscopy to measure the overlap itself. At no stage do 
we assume individual addressabihty of lattice sites, but 
rather consider global operations and the measurement 
of global quantities. 

A. Initial preparation of disorder replicas 

We consider a situation in which atoms are confined to 
a three-dimensional optical lattice, which is particularly 
deep along the 2;-direction, so that an array of planes with 
2D lattice potentials are formed in the x-y plane, with no 



tunneling between the planes. Each of these 2D planes 
of the potential can be divided into two subplanes, again 
along the z-direction, by adding a superlattice of half the 
original period which can be controUably switched on and 
off [25, 23\ ■ The resulting two subplanes then constitute 
our replicas, and our goal is that these pieces should have 
the same disorder realisation. 

In the case of disorder generated optically (using, e.g., 
speckle patterns) the different layers exhibit the same dis- 
order landscape for a suitable orientation of the lasers. 
However, for a disorder induced by a second species, the 
replicas must be produced in several steps. We start 
with a two-species mixture in the undivided plane. A 
fast ramping up of the optical lattice depth within the 
plane for the disorder species results in a Poissonian dis- 
tribution of particles in the lattice. The next step is 
then a filtering procedure removing singly-occupied sites 
and keeping only doubly-occupied and unoccupied sites 
(which are distributed randomly). Such a filtering pro- 
cedure has been achieved in recent experiments where 
atoms are combined to Feshbach molecules, and atoms 
remaining unpaired are removed from the system using 
a combination of microwave and optical fields that is 
energy-selective based on the molecular binding energy 
[55] . Finally, the superlattice is applied adiabatically pro- 
viding a double-layer structure, where exactly one atom 
of the disorder species is transferred to each of the two 
subplanes. The result is two planes with a random distri- 
bution of atoms of the "disorder species," with an exact 
copy in each subplane, i.e., a replica of the random dis- 
order pattern. This is illustrated in figure [l] 

Interactions between atoms of a probe species and the 
disorder species can now be adiabatically increased, or 
the probe species otherwise adiabatically introduced to 
the system (e.g., using spin-dependent lattices or super- 
lattices). This results in two replicas of the same system, 
which we label a and /3, with atoms in each replica evolv- 
ing independently according to their system Hamiltonian 
in the same disorder potential. We again reiterate that 
the superlattice separating the replicas should be deep, 
so that there is no tunnelling between the replicas. The 
only correlations between the two should thus be deter- 
mined by the identical disorder realisations. 



B. Measurement process 

To measure the overlap q as given by Q we need to 
compute correlation functions that compare the state of 
probe atoms in the two replicas a and (3. In each case, 
the correlation functions can be computed if we measure 
the joint probability distribution for occupation numbers 
in the two replicas (this is discussed in detail below). In 
particular, we need to determine the joint probability of 
having n particles in layer a and m particles in layer (3 
at a given lattice site, denoted by Pnm, and the probabil- 
ity of having n particles in layer a, /? at a given lattice 
site, denoted by p"'^ (both depicted in figure We 



thus first freeze the state in each replica by ramping the 
lattice suddenly to deeper values in all three directions, 
so that tunnelling no longer occurs. We can then mea- 
sure the different possible configurations of particles oc- 
cupying individual lattice sites in the layers a and P (see 
figures gHI). 

This is done in two steps: First, the layers are com- 
bined so that at each site each initial eigenstate involving 
different occupation numbers in the two replicas becomes 
a superposition of degenerate states involving the total 
number of atoms spread over motional states in a sin- 
gle well (see figure |4] section IIIB 2| below) . The final 
state is, of course, strongly dependent on the combina- 
tion process, as well as the initial state before the com- 
bination of the layers. However, we need only determine 
the fraction of lattice sites with a total of n particles, 
denoted pn, which is unchanged provided that the lat- 
tice is deep within each plane, to prevent tunnelling of 
atoms. The relative frequency with which each final con- 
figuration Pn occurs can then be measured based on the 
different spectroscopic energy shifts arising from different 
numbers of atoms in the final well, and different distri- 
butions of atoms over the motional states (see figure |5] 



section |III B 3 below). This involves weak coupling of 
atoms to an auxiliary internal state, similar to that per- 
formed in [28 . We show that the set of possible final 
states given an initial configuration can be distinguished 
spectroscopically from the set of states corresponding to 
initial configurations giving different contributions to p„, 
and thus we can reconstruct the original joint probability 
distribution pnm prior to combining the layers. The re- 
duced probabilities p^'^ can be determined from similar 
spectroscopy without first combining the layers, allowing 
the full overlap function q to be constructed. Below we 
give further details of each step in this procedure. 



1. Quantities that need to be measured 

We need to measure two quantities to determine the 
overlap q, namely [(gi)]aD, where 



(5) 



with a ^ P and L is the total number of lattice sites, cor- 
responding to the density-density correlations between 
the two layers, and [{q2)]av, [{q2)]av, where 
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(6) 



corresponding to the average density of each layer. Then 
we can express 

q ^ [{qi)U - m)]av[{q^)]av (7) 

The disorder average should be performed by repeated 
measurements of the quantum average for different real- 
isations of the disorder f60L In general, one can expect 
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FIG. 3; A sample system illustrating the required joint and 
combined system probabilities involved in the measurement 
scheme before and after combining the two layers respectively. 
Top: The individual layers a and /3 and the correspond- 
ing joint probability p„m at each site. Bottom: Combined 
layer and the corresponding combined system probability p„ 
at each site. Note that this is only a guide, because after 
combining the layers the particles are spread across multiple 
motional states (not depicted). 



that a measurement on and rij for large separation 
between sites i and j is independent of each other, and 
consequently the summation over lattice sites automati- 
cally performs a statistical quantum average [61J . In the 
following, we are interested in small filling factors such 
that lattice sites with three or more particles are rare and 
can be neglected in all three phases of the system. 

First we focus on measurement of [{q2)]av[{q2)]av, i-e., 
the second term of the overlap in (|4]) . If the layers are pre- 
pared as discussed above, the symmetry (1/L) ) = 
(l/L) ^i{n^) is preserved and thus we can measure the 



average over two replicas, [((72)]^ where |62j 
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(92 + q'2) = 



-T 



« + <)• 



(8) 



For the average density we obtain: 



(92) = ^E« 



^En(p^-l-p^), (9) 



where p'^^ — N'^'^ jL with N^'^ the number of sites with 
n particles in layers a, /3, and thus 



1 



^qi) 



(K+K + 2p^ + 2pg). 



(10) 



Next we consider the measurement of [(gi)] i.e., the 
first part of the overlap in Q . In this case the density- 
density correlations can be expressed as 



nmp„ 



(11) 
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where Pnm = Nnm/L with Nnm the total number of sites 
with n particles in layer a and m particles in layer /3, and 
thus, keeping terms up to two particles in total over the 
replicas at a given site 



(ft 



Pll + 2(pi2 +P21) +4p22- (12) 



Now consider p„ — Nn/L, with Nn the total number 
of particles from both layers, which can be measured 
by combining both layers (described in detail in sec- 
tion |III B 2 1 and subsequently applying a similar scheme 



to the one presented in PH] (described in detail in sec- 
tion III B 3 1 . Using the following equations that relate 
both p„ and p^-*^ to p„„ 



Pi 
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P12 + P21 
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Pii 



Pi 

P'l 

P2 

pI 



Pio +P12 + Pll, 
Poi +P21 + Pll, 

P20 +P21 + P22, 
P02 + P12 + P22, 



we can show that pn = ^{pi + Pi — Pi — Ps) which in 
conjunction with p^ = pi2 + P21 and P4 = P22 gives all 
necessary terms. Then, (|12|) simplifies to 



(13) 



1 3 

(91) = 2 (P? +Pi-Pi)+ + 4p4 



Finally, note that q2 as given in (10 1, can also be ex- 
pressed as (^2) = 5 Yjn^P^- 

In the following two sections we show how to mea- 
sure these conditional probabilities by combining the two 
replicas together, and performing spectroscopy where a 
single atom is coupled to a different internal state. In 



section III B 2 we describe the process of combining the 
two layers into a single one whilst in section [Tll B 3| we de- 
scribe a similar scheme to the one presented in [28^ that 
can be used to distinguish different occupation numbers. 



where Uki = 2wj^a / dz\wk{z)\^\wi{z)\'^ , with a the scat- 
tering length, the transverse trapping frequency, 
u'fe(z) the Wannier functions for band fc, and the band 
energy (fc, I S {0, 1}). The expression for Uki is a special 
case of the general 3D expression 



U = 



47r;i2 



dxdydz\wo{x)\^\wo{y)\'^\wk{z)\^\wi{z)\'^, 



where wq(x) and wo(jj) are the Wannier functions for 
the lowest motional state in the 2D replica plane {x,y), 
with the additional assumption of an isotropic and deep 
lattice potential in the replica plane characterized by the 
frequency u;±. 

States which are initially nondegenerate will evolve 
adiabatically to the corresponding eigenstate of the 
Hamiltonian in ( 14 ) , provided the barrier is lowered 



slowly with respect to the energy separation between 
the (nondegenerate) eigenstates. However, we note that 
some of the initial states are doubly degenerate due to 
the equal single particle energies in the two layers. Thus 
when the barrier is lowered these degenerate states be- 
come coupled to each other (because the single parti- 
cle states change time dependently) and the final state 
will evolve into a process-dependent superposition of the 
different corresponding eigenstates of the Hamiltonian 
in ( 14 1. However, this poses no problem for our measure- 



ment scheme since the states that arc initially degenerate 
contribute equally to the probabilities we are trying to 
measure. Moreover, for a given initial configuration it 
is not necessary to know the full final state, in fact it 
is sufficient to know only the possible corresponding fi- 
nal eigenstates of the Hamiltonian in (with the same 
number of particles as the initial configuration) that con- 
tribute to the final state. Figure|4]depicts all the different 
possible states before and after combining the layers for 
up to a maximum of four particles in total. 

In the next section we show how the different configu- 
rations for a given particle number can be measured and 
thus the total occupation numbers determined. 



2. Combining Layers 

We now describe in more detail the process of combin- 
ing the two layers into a single one and determine the 
resulting joint state for the different possible initial con- 
figurations. The combination of the layers is achieved by 
lowering the potential barrier that separates both lay- 
ers (switching off the superlattice) . Specifically, we now 
determine the final (after completely lowering the bar- 
rier) state for each of the different initial configurations 
of particles (prior to lowering the barrier). 

After the superlattice is removed and layers are com- 
bined, the system is described by the Hamiltonian 



E 

ke{0,l} 



efe^fc + -^nkin.k - 1) 



2UQinoni 



(14) 



3. Number Distribution Characterization 

We now investigate a scheme similar to the one pre- 
sented in [55] that uses density-dependent energy shifts 
to spectroscopically distinguish different particle num- 
bers at each site of the combined layer (see section HI B 2 
above) resulting in a measurement of p„. The measure- 
ment of p"'^ follows from the same spectroscopic analy- 
sis, but performed on each individual layer a, (3 before 
combining the layers. Such a scheme has already been 
implemented experimentally and was utilized to observed 
the Mott-insulating shell structure that arises due to a 
harmonic trapping potential commonly present in cold 
atomic gases experiments [28] . 

We begin with all atoms in an internal state a, and in- 
vestigate the weak coupling of particles to a second inter- 
nal state h (with at most one particle transferred). Due 
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before combination 



after combination 



tonian (within a rotating-wave approximation) 




WW 










+ 



FIG. 4: Analysis of the process of combining the two initial 
replicas, showing possible states with at most double occupa- 
tion (dots) in a given well before combination of the layers 
(left hand side), and the resulting possible states after com- 
bination of the layers (right hand side). Note that here ' + ' 
indicates a general superposition with unspecified coefficients, 
as these are dependent on the precise time dependence of sys- 
tem dynamics during combination. Knowledge of these coef- 
ficients is not necessary in the measurement scheme presented 
here. 



to the difference in interaction energy between particles 
of species a and b, and particles solely of species a there 
is an energy shift between the "initial" state (all particles 
of species a) and "final" state (one particle of species b 
and all others particles of species a) . Thus different ini- 
tial particle numbers can be distinguished by the energy 
shift (corresponding to the Raman detuning at which 
the transfer is resonant), provided the shift is different 
for the different particle numbers. Since particles in the 
combined layer can be spread across multiple motional 
states (see section III B 2 above) , different energy shifts 



occur for different configurations of the same number of 
particles. Whilst we find that the energy shifts of some 
of the configurations (of the same number of particles) 
are identical or very similar, there are other configura- 
tions having substantially different energy shifts. Thus 
we have to check that the energy shifts for all configura- 
tions of a given number of particles can be distinguished 
from energy shifts of configurations corresponding to a 
different total number of particles. To this end we now 
present a detailed calculation of the energy shifts arising 
from all the different initial configurations of particles 
and show in which cases the total number of particles 
can be distinguished. 



Let aj, and b], denote the creation operators for par- 
ticles of species a and b, respectively, in band k where 
k G {0,1}. Then the Raman coupling between the two 
internal states a and b is described by an effective Hamil- 



Hnc = ^{albo + a\bi^ 



^R.c.)~A{t){blbo + b\bi), (15) 



where A{t) is the Raman detuning, ri(i) the effective 
two-photon Rabi frequency, and we have assumed that 
ri(t) ^ A' where A' are the detunings of the lasers from 
the atomic excited state. In what follows we assume that 
J Q,{t)dt ^ TT, i.e., weak coupling to internal state b and 
that the two lasers creating the Raman transition are 
running waves with equal wave vectors. Note that pro- 
cesses in which the internal state and the band would 
change do not occur provided the Raman detuning and 
effective two-photon Rabi frequency are smaller than the 
band separation. 

The onsite interaction energy for atoms of species a and 
6 in the lowest two bands is described by the following 
Hamiltonian 



E 

ie{a,6} 



TTll 

"^00 

2 



+2C/-n^nl 



Tjii 

"oK-i) + ^"lK-i) 



-n°no -I- {a\aib\bQ + a\a^)b\bi + H.c.) 



where 



2uj±aij / dz\wk{z)\'^\wi{z)\'' 



(16) 



(17) 



•^00 



with Qij the scattering length between two particles in 
internal states i and j G {a, fe}), lj±_ the trans- 

verse trapping frequency, and Wkiz) the Wannier func- 
tions corresponding to band k {k,l £ {0,1}). Note 
we have also introduced the particle number operators 
"fc — o.\ak and n^. — b\bk- We further make the as- 
sumption that the lattice is very deep so that we can ap- 
proximate all the Wannier functions by simple harmonic 
oscillator eigenfunctions. This assumption of a deep lat 
tice gives the following simplified expressions 2U^\ 

and U'A = (3/4)C/o^o- 

Using ( |16[ ) we now explicitly compute the energy shift 
AE — Ef ~ Ei associated with the processes \i) |/), 
where \i) is the initial state corresponding to all particles 
of species a and having energy Ei, and |/) is the final 
state corresponding to one particle transferred to internal 
state b and having energy Ef. Clearly, AE depends on 
how the particles are initially distributed in the lower and 
higher bands after combining the two layers (see the right 
hand side of figure |4]) . There are two different cases to 
be considered; either all particles are initially in the same 
band (this is the situation in |28j ) or particles are initially 
in both bands. In the latter case it is then possible that 
particles in either band are transferred to internal state 
b and since these different possible resulting states are 
coupled by the interaction Hamiltonian in ( 16 1 the final 
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resulting state is a superposition of these, corresponding 
to an eigenstate of the Hamiltonian in (16 1. 



initial 



->• final 



aa 
00 



For example, let us consider a total initial number 
of two particles as depicted in figure [Sj The different 
initial and final states are most conveniently expressed 
in the occupation number basis \Nq , N^; Nq, N^) with 
j^a,b j-^o.fe^ ^Yie number of particles of species a, b in the 
lower (higher) band at a given lattice site. First con- 
sider the case where both particles are initially in the 
lower or higher bands corresponding to the initial states 
|2,0;0,0) or |0,2;0,0) (see top and bottom of figure [s]) . 
The corresponding final states are given by |1,0;1,0) 
and |0,1;0, 1), and the corresponding energy shifts are 
AE/Ug^ = e- 1 and AE/Ug^ = f (e- 1), where we have 
defined e = Uqq/Uqq. Next we consider the case where 
there is a particle in each band initially corresponding 
to the initial state |1,1;0,0) (see middle of figure |5]). 
Then the corresponding final states are the eigenstates 
of the Hamiltonian in ( 16 1 in the subspace correspond- 



ing to one particle of species a and one particle of species 
& given by |1, 1)± = (1/^2) (|0, 1; 1, 0) ± |1, 0; 0, 1)), with 
eigenenergies E+/Uqq = e and E-/Uqq = 0. The associ- 
ated energy shifts for these two possible final states are 
AE/UgS = e - 1 and AE/U§g = -1 respectively. Note 
that the Raman transition only couples the state |1, 1) + 
to the initial state |1, 1; 0, 0) with a matrix coupling el- 



while for the 



ement fi/i = •s/Z, where ilj^ 

state |1,1)- we find fi/i = 0. Figure [s] summarizes the 
different possible configurations for a total of two parti- 
cles together with the corresponding energy shift AE and 
matrix coupling element f2/i. We note that the energy 
shifts, AE, for the different final states to which the ini- 
tial states couple (i.e., ^0) are very similar, this is 
of advantage since our scheme needs only to distinguish 
energy shifts for different total particle numbers. 

Similarly, we can determine the energy shifts and cou- 
pling matrix elements for other numbers of particles. If 
all particles are initially in the lower or higher state corre- 
sponding to the initial states |7V^, 0; 0, 0) or |0, iVf ; 0, 0), 
the corresponding final states are given by jiVp — 1, 0; 1,0) 
and lOjA'^f — 1;0,1). The corresponding energy shifts 
are then given by AE/Ugg = {Ng - l)(e - 1) and 
AE/Ugg = |(iVf - l)(e - 1), and the matrix coupling 
element in each case is ^Ifi = s/Ng and {Ifi — ^/Nf 
respectively. For other initial configurations with m par- 
ticles in the lower band and n particles in the higher band 
the corresponding final states are found by diagonalizing 
the Hamiltonian in (16) in the appropriate particle sub- 



space. This always results in a pair of eigenstates denoted 
\m,n)± with eigenenergy and using these expres- 

sions the energy shift AE and coupling matrix elements 
n fi can be directly determined. For a total initial num- 
ber of three and four particles the expression s for |m, n)± 

(m,n 



E. 



AE and f2 fi are given in appendix 



Al 



and 



A2 



respectively, for all the different initial configurations. In 
table |T] we list the results for up to a maximum num- 
ber of four particles; the first column lists the different 



flj; AE/U, 

\/2 (e-1) 
V2 (e-1) 



V2 



(e-1) 



FIG. 5: Illustration of the number distribution characteri- 
zation within the two layers for states with two particles in 
the combined layer. On the left are all possible initial config- 
urations of two particles of species a in the combined lattice 
(dark-coloured circles) and the corresponding final states with 
one particle each of species a and b (light-coloured circles). On 
the right we list the respective coupling matrix elements, fi , 
and energy shifts, AE/Uqq, for each of the configurations on 
the left. 



initial configurations and corresponding final states, the 
second column lists the corresponding matrix coupling 
elements fi/i and the third column lists the associated 
energy shifts AE. 

It is possible to distinguish different occupation num- 
bers for specific ranges of e if the energy shifts for different 
total numbers of particles are distinct. The desired range 
of e is then between points where any energy shifts for 
different total number of particles would coincide (e.g., 
e = 1) For a specific example we consider the value 
e — 0.98 (corresponding to scattering length values in 
|28j ) and give numerical values for the matrix coupling 
elements Ufi and the energy shifts AE in columns four 
and five, respectively, in table [Tj Then, to a high accu- 
racy, the determination of p„ is obtained by applying the 
coupling between the different internal states a and b for 
a short time t„ = At / ^Jn at all resonant frequencies in 
table |l] with total number of particles n (ignoring the res- 
onances with small and vanishing couplings f2/i). The 
total number of transferred particles is then proportional 
to p„ and efficiently avoids problems arising from non- 
adiabatic combination of the layers. 



IV. OVERLAP FUNCTION IN THE 
DISORDERED BOSE-HUBBARD MODEL 

In this section we determine the overlap function q 
[see ([4])] in the disorder Bose-Hubbard model, both an- 
alytically and numerically. We focus on a bimodal dis- 
order potential = ±A with the probability distribu- 
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TABLE I: Resonant energies and coupling strengths for the 
transfer of single atoms in the combined layer to a different 
internal state. Column one lists the different initial configu- 
rations and the corresponding final states up to a maximum 
number of four particles. In the second and third column we 
list the matrix coupling elements Qfi and energy shifts AE, 
respectively, corresponding to the different initial configura- 
tions in the first column. The states |m,n)± are listed in the 
text and appendices together with the corresponding eigenen- 
ergies E^'"'' and coupling matrix elements r;j!"'"'- The last 
two columns give numerical values of Sl/i and AE for the 
specific value e — 0.98. 
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tion P(A) = 1 — P(— A) — p, which is the appropri- 
ate description for the chosen disorder implementation 
(see discussion in earlier sections). In the following, we 
present the modifications due to the presence of a weak 
(A < U/2) bimodal disorder potential, and the predic- 
tion for the overlap function within the different phases 
(see figure [2| . It follows from this constraint that we do 
not enter the compressible disordered phase, and hence 
the system will always have a unique ground state. Thus, 
the ground states in two exact copies a, /3 of the system 
(i.e. same disorder distribution) are the same, and conse- 
quently (nf ) = (nf ). We therefore expect a Sherrington- 
Kirkpatrick-type behaviour for the overlap [64j . i.e. the 
overlap between any two copies will always be the same, 
and thus the overlap function q can in fact be calculated 
using 

While the above is a self-evident statement for a 
replica-symmetric system such as the one we examine 
here, it may no longer be true for systems thought to 



exhibit replica-symmetry breaking (e.g. a classical spin 
glass type model [58j , or quantum systems with extended 
interactions [Si]). For the latter type of system, measure- 
ment of the overlap function using physical replicas (see 
section III) is thought to yield different values for the 
overlap depending on a and /?, the defining signature of 
replica symmetry-breaking. 

In section |IV A| we analytically calculate the overlap 
function, for the various phases of the disorder Bose- 
Hubbard model ([ij, whilst in section IV B we present 
results for the overlap function from numerical simula- 
tion of the one-dimensional system. 



A. Analytical determination of the overlap function 

Starting with the limit of vanishing hopping J = 0, ([iJ 
reduces to an onsite Hamiltonian and the ground state is 
given by jfi) — 11^151)^ with the lowest energy state 
within each site i. This lowest energy states takes the 
form |f2)i = \ni) with \ni) being a state with fixed particle 
number rii at site i, with Ui subject to the constraint 



J7 (n^ - 1) < ^ - Ai < Urii 



(18) 



For a weak bimodal disorder potential (A < U/2), we 
have to distinguish two different cases. 

First, for a chemical potential within the range U^uq — 
1) -I- A < /i < Uno — A the above constraint in (18 1 
is fulfilled independently of the site parameter i for the 
integer value uq, i.e., at each site the ground state is 
characterized by the integer particle density [(ni)]o„ — uq 
and we obtain a Mott-insulator (MI) phase (see figure |2|. 
Furthermore, the overlap parameter vanishes identically 
in this limit, i.e., q = for a Mott-insulator at J = 0. 

Second, for chemical potentials in the range Uhq — A < 
fi < Uno + A, the particle number at sites with differ- 
ent disorder potential A,; differs by a single particle, i.e., 
at each site with disorder potential A^ = A the particle 
number is determined by no, while at sites with disorder 
potential A^ = —A the lowest energy state is character- 
ized by riQ -l- 1 particles. Therefore, the averaged particle 
density and the overlap parameter take the form 



[{ni)]av = n„ + (1 -p) 



q = p{l -_p). 



(19) 



This situation corresponds to the disorder-induced Bose- 
glass (BG) phase (see figure |2]) and the ground state is 
characterized by nonhomogeneous filling. Similarly to 
the case of the MI phase, the particle number is fixed 
over a finite range of the chemical potential and thus 
the BG is incompressible. In contrast to the MI phase, 
the BG phase is characterized by a non- vanishing overlap 
parameter dependent on p. 

Next, we focus on small hopping, J <C J7, A, /x, and de- 
termine the overlap function using perturbation theory. 
We seek a unitary transformation of the form U — e~*'^ 
which transforms the total Hamiltonian in ([l| to an ef- 
fective Hamiltonian having the same eigenspectrum but 
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acting only in the space of the unperturbed eigenstates. 
Using this unitary transformation we can directly deter- 
mine the corrections due to hopping in perturbation the- 
ory by computing (ri|ni|ri) where fij = U^riiU. Details 
are given in appendix [B] 

We start with the Mott-insulator and consider the 
leading order (non-vanishing) correction to the overlap 
function which occurs at fourth order due to a site 
(disorder) independent density, (il|ni|ri) = tiq (using 
= l^o) in the MI). To lowest order in perturbation 
theory the overlap is then given by 



where 



with 



^2 = -2^Si^ [Si,ni\], 



Si = -i'^Cjk\i^){(f>jk \ +H.C. 



(20) 



(21) 



(22) 



and \(j).jk) = \no + l)j|no - l)fcn;^y |rio)/. Using these 
expressions it is straightforward to compute the expecta- 
tion value of the density correction 



(ia) = no{no + 1) 



Uii)) 



(23) 

where E^ — Eij = —U +Aj— A^. It is now straightforward 
to calculate the disorder averages and the overlap is given 

by 

q = 64(no(no + l))Ml " P)<z + (24) 

where z is the number of nearest neighbours. 

Next we turn to the Bose-glass phase and consider the 
leading order (non-vanishing) correction to the overlap 
function which occurs here at second order. The expres- 
sion for the overlap in lowest-order perturbation theory 
is then given by 

q=qO + J^ {[{n^){A2)]av - [{n^)]av[{A2)]av) , (25) 

where qo = p{l — p) (the overlap at zero hopping), A2 
and 5*1 are as given above for the MI phase but with 

\^) = n [^'1^0 + + (1 - di)\no)i] (26) 



and 

\(t>jk) = [dj\no + 2)j + {1- dj)\nQ + [dk\no)k 
+(1 - dk)\no - l)k]^i^jk[di\no + 1)/ + (1 - di)\no)i], 

where dj k^i (1 — dj^k,i) is zero (one) for Aj j, ; = A 
{Aj^k.i = ~A). Again, using these expressions it is 



straightforward to compute the expectation value of the 
density correction and for no > 1 we obtain 
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{Eo-E,,Y iEo-E,,y 



(27) 



where 



and 



Eo - E,, = -C/(l - d, + dj) + Aj - A, (28) 



= {d,V^^^ (29) 
+(1 - rf,)Vno + l)(djVno + l + (l - dj)y^). 

Once again a straightforward calculation of the disorder 
averages gives the following expression (tiq > 1) for the 
overlap: 

n ^f^ ^2 (^0 + l^UjU - 2A) + 
q-p{l~p) [l^J 2A^{U-Ar )■ ^^^^ 



An analogous calculation for no = gives 



(31) 



which corresponds to the U limit of ( 30 1 



When the effects of hopping dominate, the bosonic 
atoms are in the superfluid phase (see figure 2]). Then, 
the influence of weak disorder on the superfluid phase 
can be studied within a generalized mean-field approach. 
In analogy to the Gross-Pitaevskii formalism, we replace 
the bosonic operators bi in the Hamiltonian (jlj by a lo- 
cal complex field V'i , and minimize the corresponding free 
energy functional H(ipi). The effect of disorder is to pro- 
duce fluctuations in the local density. Hence we make the 
ansatz V'l = e^^V'T'O + Sui where no is the homogeneous 
density in the absence of disorder and Srii ex A^ is the 
disorder induced density fluctuation. Introducing the no- 
tation Sn — [Sn.i\av as the disorder-averaged density fluc- 
tuation, the particle density reduces to [(ni)]ot, — nQ+6n. 
We substitute this ansatz into H{ipi), expand tpi for small 
fluctuations, Srii <C no (see appendix C| and minimize 
the resulting expression in ( |C1| at flxed particle number, 
which is equivalent to requiring [^ni]c„, = 0. Note that 
for large dimension and thus a large number of nearest 
neighbours we may replace 



Uii)) 



zSn. 



(32) 



Using this relation and self-consistently solving for the 
density fluctuations Sui (see appendix [C]) gives the over- 
lap function 



q = 



( 4710 A 



\2UnQ + zJ 
where no = (/i -I- Jz)/U + ^. 



P(l -P), 



(33) 
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In summary, we find that in the limit of small hopping 
J /U <C 1, the overlap signals a sharp crossover between 
the Mott-insulator with g « and the Bose-glass with 
q Ki p{\ — g), see figure |2] In turn, the superfiuid phase 
is characterized by off-diagonal long range order giving 
rise to coherence peaks in a time of flight picture. Con- 
sequently, the overlap and the coherence peaks in time of 
flight allow for a clear identification of the phases in the 
disordered Bose-Hubbard model. However, we would like 
to point out that the overlap function behaves smoothly 
across the phase transition and consequently is not suit- 
able for the precise determination of the exact location 
of the phase transition. 



B. Numerical results 

We have simulated the disordered Bose-Hubbard 
model in ID using a number-conserving algorithm 
based on time-dependent matrix product states [551 167j . 
In the matrix product state representations, we reduce 
the Hilbert space by retaining only the x most significant 
components in a Schmidt-decomposition of the state at 
each possible bipartite splitting of the system. States in 
ID can be efficiently represented in this way, with rel- 
atively small X giving essentially unity overlap between 
the represented state and the actual state. The simula- 
tions both serve as a check for the perturbation theory, 
and allow computation of the overlap function beyond 
the perturbation limit from the microscopic model. 

Because we consider the regime 2 A < [/, the sys- 
tem remains incompressible and does not enter the dis- 
ordered glass phase. Therefore we do not expect replica- 
symmetry breaking to occur. For our simulations, the 
presence of replica-symmetry is equivalent to convergence 
to a unique ground-state for any given disorder indepen- 
dent of initial conditions. The overlap function can then 
be computed as 



9= J^^{n^)'^ -n^ 



(34) 



where n is the density of the system. Computation of 
the overlap function can be simplified in this way as 
L~^^-(ni) self-averages to [{ni)]av for sufficiently large 
systems, which in turn self-averages to n. 

The ground states used in p4| are computed by 
imaginary-time evolution, c.f. |67j . of the ID Bose- 
Hubbard-model given by ([T]) with a randomly-chosen bi- 
modal distribution, and different impurity strengths 2A. 
The lattice size is L = 64, representative of current ex- 
perimental capabilities, and we compute the overlap as 
a function of A for several different values of U / J and 
the fraction of impurity-occupied sites, p. In figure [6] 
we display results for U/J= 10 and p = 0.5. We have 
performed convergence tests in the time step, the dura- 
tion of imaginary-time evolution, the number of disorder 
realizations, and the number of retained states %. We 



find that the overlap function converged even for sur- 
prisingly small values of x ~ 20 states in the large U / J 
limit (Mott-insulator or Bose-glass regime). Regarding 
the number of disorder realizations for each value of A, 
we find that whilst 10 random realizations of bimodal 
disorder are not sufficient to narrow the error down, 20 
realizations suffice. Much of the observed fluctuations in 
the value of the overlap function for a given A stem from 
the fact that we allow the number of impurities in each 
realization to fluctuate around the mean Lp. Conversely, 
in our simulation data we observe different disorder real- 
izations with the same number of impurities having over- 
laps closer to each other than the error bars is figure [6] 
would suggest. 

As shown in figure [6] we generally found very good 
agreement between numerical simulations and pertur- 
bation theory in its region of validity. For the Mott- 
insulator regime (commensurate filling factor) this is the 
regime where 2A <C U . We observe the breakdown of 
our second-order perturbation theory as A approaches 
C//2 — 5 J, which is caused by energy-degeneracy of ad- 
jacent impurity- and non-impurity occupied sites when 
A — U/2. In the Bose-glass regime, breakdown of per- 
turbation theory sets in for low A, as degeneracy occurs 
there for A = 0. 



0.15 




0.05 



FIG. 6: Overlap function g as a function of impurity strength 
A in the Mott-insulating phase at filling factor uq — 1 (black 
lines) and in the Bose-glass phase for tiq = Q (grey lines). In 
both cases, the impurity probability is p = 0.5 and U / J = 10, 
for L = 64. Simulation results (solid lines) are shown versus 
results from perturbation theory (dashed lines). Much of the 
variation in the simulated overlap function is due to the fact 
that we allow the number of impurities to fluctuate around the 
value given by Lp in each particular disorder realization. Note 
that the perturbation result in the Mott-insulato r di verges as 
A approaches U due to energy degeneracy [c.f. (24l]. In the 
Bose-glass phase, second-order perturbation theory yields a 
divergence as A approaches 0, again due to energy degeneracy 
[c.f ((3l|]. 
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V. CONCLUSION 

We have shown how the Bose-glass phase can be iden- 
tified by measuring the overlap function characterizing 
the correlations between disorder replicas, i.e., systems 
having identical disorder landscape. We have described 
a procedure to create disorder replicas using cold atomic 
gases in an optical lattice, focusing on the case of dis- 
order induced by a second atomic species. A scheme to 
measure the overlap has been presented, which involves 
the characterization of the occupation numbers in both 
the individual replicas and the joint system, where both 
replicas have been combined. Specifically, we have shown 
that after combining the replicas together, particles can 
be distributed across multiple motional states and we ex- 
plained in detail how to determine the occupation num- 
ber distribution in this situation. Using perturbation 
theory we have calculated the overlap for weak hopping 
and have shown the different behaviours exhibited in the 
Mott-insulator and Bose glass phases; these results are in 
good agreement with the results from a one-dimensional 
numerical simulation. In the opposite limit of very large 
hopping we also obtain analytical results for the overlap 
within a mean- field-theory treatment. 

Applying our proposed measurement scheme for the 
overlap to other types of disordered models [551 Ei] would 
provide interesting insights to characteristic properties 
of disorder-induced quantum phases, such as the pos- 
sible appearance of broken ergodicity and different de- 
generate ground states. Of particular interest would be 
the quantum spin glass where the overlap corresponds to 
the Edwards-Anderson order parameter which has been 
conjectured to exhibit replica symmetry breaking |58j . 
i.e., the results for the overlap depend on the particu- 
lar replica. Our method could also give novel insights 
into physics of systems that possess multiple metastable 
states, such as ultracold dipolar gases US] or quantum 
emulsions |69j ; in particular statistics of overlaps of such 
metastable states could be measured. In addition to equi- 
librium properties considered here, the dynamics of the 
overlap might also posses an interesting behaviour, and 
could be measured in a similar setup. In one-dimension 
the time evolution of the overlap could be determined 
numerically [67^ . 
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APPENDIX A: ENERGY SHIFTS AND 
COUPLINGS FOR THREE AND FOUR 
PARTICLES 

In the following two sections we determine expressions 
for the eigenstates and eigenenergies of ( 16 1 , i.e., \m, n)± 



and E^^^ and associated energy shifts and couplings 
for all possible initial configurations of three and four 
particles. Note that the case where all particles are in the 
lowest or highest energy band has already been explained 
in section ITlI B 31 and will not be discussed here further. 



1. Three Particles 

We consider the two different configurations described 
by the initial states i) |2, 1;0,0) and ii) |1,2;0,0). 
For case i) the two eigenstates are given by 

|2,1)+ = -^(|2,0;0,l)+^/2|l,l;l,0)), (Al) 

|2,1)_ = -^(-V2|2,0;0,l) + |l,l;l,0)) (A2) 

with eigenenergies E\_ = 2e -t- 1 and E_ ' = e/2 -t- 1 
respectively. We see that only the state |2, 1)+ couples to 
the initial state |2, 1; 0, 0) , with a matrix coupling element 
of flfi = \/3, whilst for the state |2, 1)_ we find flfi = 0. 
The energy shifts for the two different possible processes 
are given by AE/Ugg = 2(e - 1) and AE/Ug^ = f - 2. 
For case ii) the two eigenstates are given by 



|1,2). 



where 



(7£'"|0,2;1,0) 



|i,i;0,i; 



(1,2) 

7± 



-1 - e ± \/33e2 + 2e + 1 



4^/2e 



(A3) 



(A4) 



and with eigenenergies 



E. 



(1,2) _ 



e-f l±^v/33e2-H2e-hl4- ^(e- 1). 



(A5) 



In this case both states |1, 2)-|- have nonzero coupling to 
the state |1, 2; 0, 0) given by 

flf, = Ty^-') EE (l,2|±i7«c|l,2;0,0) (A6) 
= , \ , (7i^'^^ + V2). 

The energy shifts for the two different possible processes 
are given by AE/Ugg = E^^''^'' - 3. 
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Four Particles 



We consider the three different configurations de- 
scribed by the initial states i) |3, 1; 0, 0), ii) |2, 2; 0, 0) and 
iii) |1,3;0,0). 

For case i) the eigenstates are given by 



|3,1> + 
13,1)- 



1 

1 

7i 



(|3,0;0,1> + V3|2,1;1,0)), (A7) 
(-\/3|3,0;0,l> + |2,l;l,0)), (A8) 



with eigenenergies E^^'^^ — 3(e + 1) and iJ^'^'^^ — g + 3^ 
respectively. We see that only the state |3, 1)+ couples to 
the initial state |3, 1; 0, 0), with a matrix coupling element 
of ri/i = -x/i, whilst for the state |3, 1)_ we find fi/i = 0. 
The energy shifts for the two different possible processes 
are given by AE/U^o = 3(e - 1) and AE/Ug^ = e - 3. 
For case ii) the eigenstates are given by 



|2,2)d 



where 



(7i''|l,2;l,0) + |2,l;0,l) 



(A9) 



and with eigenenergies 



E. 



(2.2) _ 



3 + 2e ± \/65e2 - 2e + 1 - ^(1 + e). 



(All) 



In this case, generally both states |2,2)± have nonzero 
coupling to the state |2, 2; 0, 0) given by 

= ?7±^'^ = (2,2|±iJ:,.c|2,2;0,0) (A12) 



The energy shifts for the two different possible processes 
are given by AE/Ug(^ = e'^''^^ - f . 

For case iii) the eigenstates are given by 

|1,3)± = , \ (7£'''|0,3;1,0) + |1,2;0,1)) , 

(A13) 



where 



and with eigenenergies 



E. 



(1,3) 



3 + 2e±\/l3e2 + 2e+l- 7(2 + 6). 



(A15) 



In this case, generally both states |1,3)± have nonzero 
coupling to the state |1, 3; 0, 0) given by 



/. = r4i^3' = (l,3|±i7,c|l,3;0,0) (A16) 
(7r^ + V3). 



The energy shifts for the two different possible processes 
are given by AE/Ug^ = E^^'^^ - f . 



APPENDIX B: PERTURBATION THEORY 

In this section we derive the corrections to the over- 
lap function due to small hopping in perturbation theory 
following closely the method presented in [70]. We start 
with the zero hopping Hamiltonian Hq with ground state 
and eigenenergy Eq. We treat the hopping Hamil- 
tonian, JHi with Hi = — which has eigen- 
states denoted by and eigenenergies denoted by Eij 
in perturbation theory. We proceed by determining a 
unitary transformation that transforms the total Hamil- 
tonian H = Hq + J Hi into an effective Hamiltonian H' 
with the identical eigenspectrum but acting only in the 
unperturbed Hilbert space of the state |f2). The required 



unitary transformation may be written as U 



and 



can be determined consistently to any order by writing 



S ^ JSi + J'S2 



We assume that S has only 



non-diagonal matrix elements which connect the unper- 
turbed ground state \Q) and the excited states Us- 
ing the unitary transformation we can directly calculate 
the correction to the particle density due to the hopping 
by computing (ri|ni|ri) where hi = U^riiU. We show in 
the following sections that it is sufficient to determine S 
to first order. Moreover, it is straightforward to show 
that to first order the matrix elements of S, i.e. ^i, are 
given by 



En — Ei. 



(Bl) 



1. Mott-insulator 

In this section we derive the leading order correction 
to the overlap function for small hopping in perturbation 
theory for the MI phase. We begin by expanding S to 
fourth order and calculate the average corrected density 
to be 



4 

E 



(B2) 
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where and 
Ai = (B3) 
A2 = [S2,n,]-]^[S,,[S,,n,]], (B4) 

(B5) 

Aa = [5'4,ni] - i[S'2, [S'2,ni]] + ^[S*!, [S*!, [S*!, [S'i,ni]]]] 



i([5i,[53,n.]] + [^3,[^i,ri,]]). 



(B6) 



First, note that since we have chosen S to have no diag- 
onal matrix elements ([5'„,ni]) = to all orders. 

We now calculate the first part of the overlap function 
as given in ([s]) and find 

[{n^)\v = [{n^)\v + J\n^){A2)]av + j\ni){A^)]av 

+J%{n,){AA) + {A2)X, + 0{j''). (B7) 

Similarly, for the second term of the overlap in (|3| we 
obtain the following expression 

[(".)]L = [{n^)]L + J\n,)U[{A2)U m 

av 

+ J^{[{n^)]av[{AA)]av 
+ [(A2)]L) + 0(J^). 

Using {n.i) — uq we obtain the following expression for 
the overlap 

= J*([(A2)^].„-[(^2)]L)+0(J')- (B9) 

Note that because (n^) = uq at all lattice sites (indepen- 
dent of the disorder) the second- and third-order contri- 
butions cancel. 

To determine the required commutators [Si, [Si, n^]] 
we write Si and rii as projectors 



Si = -i^cjfc|f})((/)jfc| +H.C., 

rij = no\n){n\+nQ ^ \<Pjk) {(f'jkl 

0«> 



(BIO) 



where 



•jk) 



(B12) 



Eq — Ejk 
\n)=l[\no},, (B13) 



\^jk) = \no + l)j\na - 1),. [] \no)i. (B14) 

Using these expressions it is straightforward to compute 
the expectation value of the density correction 



{A2) ^ noino + 1) (jW^ 
Uii)) ° 

where Eq - Eij = -U + ~ A 



(B15) 



2. Bose-glass 

In this section we derive the leading order correction 
to the overlap function for small hopping in perturbation 
theory for the BG phase. Again we consider the corrected 
density (fii) as given in (B2) but only up to second order 
in J. Again we have {[Sm rij]) = because we have 
chosen S to have no diagonal matrix elements. Thus the 
overlap to second order is given by 

q = [{nrf]av - [(ni)]L 

= qo-J^ ([(«,) {A2)U " [{n^)]av[{A2)]av) + 0{J^), 

(B16) 

where go = [{ni)'^]av — [( n-.,)]^^ is the overlap at zero hop- 
ping derived in section IV A Since (n^) is now site de- 
pendent the second-order contribution does not vanish. 

To determine the required commutators [5*1, [5*1, n^]] 
we again write Si and rii as projectors 



Si = -«^CjA.|f^)(0jfc| +II.( 



(B17) 



= {no + di)\n){n\ + {no + di) ^ |'/>jfc)(<?!>jfc| 

(jk)^t 

+ {d,{na -I- 2) 4- (1 - d,)(no + 1)) ^ |0.y)(<?!>iil 

OX*)) 

+ (d,no + d,){no ~ 1)) \'l>3^)i^J^l (B18) 



where 



Cjk 



Eq - Ejk 



(BU) with 



P) =\{[di\no)i + {l~di)\no + l)i\ 



(B19) 



(B20) 



and 



\(j)]k) = [dj\n(i + 2)j + {I ~ dj)\no + l)j][dk\no)k 
+ (l-4)|no-l)fc] 

xn,^,fcK|no + l)i + di)\no)i] (B21) 
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where dj,fc,/ (1 — dj^k.j) is zero (one) for Aj^k,i = A 
{Aj^k.j = —A). Again, using these expressions it is 
straightforward to compute the expectation value of the 
density correction and for no > 1 we obtain 



Next, we minimize (CI I at fixed particle number which 

= 0. 



is equivalent to requiring [Sui 
consistently solving gives 



Using ( 32 1 and 



{A2 



where 




7,, 



(Eo - E,,Y 



(B22) 



Srii = 



A{2p- 1) - A 
U + — 

^ ^ 2no 



(C2) 



jij = (djV^io~+2 + (1 - di)y/no + l) (d^Vf^o + 1 

+ (1 - (B23) 

and Eq - E.ij = -U[di + (1 - dj)] + A^ - A^. 



Finally, then the overlap can be expressed in terms of 
these fluctuations as follows 



APPENDIX C: MEAN FIELD THEORY 



Wc replace the annihilation operators bi in the Hamil- 
tonian in ([l]) by the mean field '(/'i = 6**^-^/^0 + ^"^^i ^^^d 
expand the fiuctuations Srii <C up to second order to 
obtain 



H{5n,) 



{5n] 



' -. — SuiSrij — - 
4nof- ' ' 2 

i 

-^{[1 + zJ - AC)8n, + 0{8n\). (CI) 



= {^'AXav - [(5? 



2 

av ■ 



(C3) 



Substituting ( C2 1 into ( C3 1 and performing the dis 



order average for a bimodal distribution gives the result 
in ( 33 1 for the overlap function q. 
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